Multi-mode description of an interacting Bose-Einstein condensate 
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We study the equilibrium dynamics of a weakly interacting 
Bose-Einstein condensate trapped in a box. In our approach 
we use a semiclassical approximation similar to the descrip- 
tion of a multi-mode laser. In dynamical equations derived 
from a full TV-body quantum Hamiltonian we substitute all 
creation (and annihilation) operators (of a particle in a given 
box state) by appropriate c-number amplitudes. The set of 
nonUnear equations obtained in this way is solved numeri- 
cally. We show that on the time scale of a few miliseconds 
the system exhibits relaxation - reaches an equilibrium with 
populations of different eigenstates fluctuating around their 
mean values. 



Bose-Einstein condensation of an ideal gas is a fa- 
mous example of a phase transition which relies purely 
on quantum statistics. The recent experimental achieve- 
ment of Bose-Einstein condensation in a dilute gas of 
trapped alkali atoms [|l|-Q triggered a revival of inter- 
est in various properties of the condensate. Some old 
but fundamental problems of statistical description of the 
Bose-Einstein condensate have been extensively studied. 
In particular, the issue of microcanonical and canonical 
fluctuations of a nonintcracting condensate was solved 
successfully |5|-|lT| . On the contrary, a closed system de- 
scription of fluctuations of the weakly interacting con- 
densate seems to be still a controversial problem at least 
at temperatures close to the critical one fl^-[l5|. The 
standard mean field description based on the Bogolubov 
approach fails |^6| at the critical temperature, as then 
the order parameter vanishes. On the other hand many 
excited states become highly occupied particularly in the 
case of realistic systems consisting of a finite number of 
atoms. The complexity of the problem grows further 
when one is interested in an estimation of the charac- 
teristic time scale of the condensate fluctuations [|l^,|l8| . 
The calculation of the two-time correlation function re- 
quires obviously a dynamical approach. 

In the following, we present an approach that allows 
the study of the time evolution of an interacting Bose 
gas. The problem is very difficult and existing approaches 
lead to numerical algorithms which are hard to imple- 
ment. The first one [g^ 20 1 invokes quantum kinetic the- 
ory and leads to a quantum Boltzmann master equation. 
The second one [^,^ is based on the separation of the 
"classical" mean field describing the condensate from its 
fluctuations. The fluctuations are quantal in nature and 
the dynamics couples averages of both their normal and 



anomalous parts to the condensate mean field. None of 
the two mentioned methods seems to be capable of han- 
dling a realistic case of a large number of particles at a 
relatively high temperature. 

The first problem one encounters in describing the in- 
teracting condensate is the real N-body nature of the sys- 
tem. The link between this problem and the traditional 
"giant matter wave" description is not straightforward. 
The only exception is the case of a gas trapped in a box 
with periodic boundary conditions. Here the symmetry 
of the potential uniquely enforces the form of the eigen- 
states of the one-particle density matrix regardless the 
interaction. Due to the symmetry, the condensate can 
be associated with the ground state of the trap. There- 
fore, to avoid all ambiguities related to the identification 
of the condensate we are going to study here the system 
in a box with periodic boundary conditions . 

We want to present our approach rooted in the theory 
of multi-mode lasers. Let us start with a general for- 
mulation of the N-body problem. The second-quantized 
Hamiltonian for the atomic system confined in a box with 
periodic boundary conditions and interacting via pair- 
wise contact potential may be written in the following 
form: 
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where $ is an atomic field operator, V = is a. volume of 
the system (L being a size of the box) and g — 4?Tr^ char- 
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acterizes the atom-atom interactions in the low-energy, 
s-wave approximation (a being the scattering length and 
m - the mass of the atom). The field <& is expanded in 
natural modes of the system - the plane waves: 
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where at are bosonic annihilation operators, and k ~ 
with rii = 0,±1,±2,... (i = x,y,z). With this 
substitution the Hamiltonian assumes its final form: 
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where £ = J^(^)'^- After elimination of a fast time 
dependence with the substitution at = exp(— i^n^t)ak, 
the Heisenberg equations of motion for the operators ak 
acquire the following form: 
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Solving the nonlinear operator equations is diffi- 
cult. The complexity of the problem obviously requires 
some approximation. A semiclassical approximation is 
particularly well suited for the description of a finite sys- 
tem at temperatures below the critical one, except the 
region close to the absolute zero. 

The semiclassical approximation consists in replacing 
all operators ak by c-number complex amplitudes (we 
are not going to introduce a separate notation for corre- 
sponding complex fields). At very low temperatures only 
the lowest lying states are macroscopically occupied and 
quantum fluctuations in excited states become important 
(see [0 for comparison). Therefore, the relevant range of 
temperatures for the applicability of our model excludes 
temperatures close to zero. From the viewpoint of the 
Bogolubov method such an approach is legitimate 
as indeed many modes are macroscopically populated, i.e. 
their occupation is greater than quantum fluctuations. 

The semiclassical approximation leads to nonlinear dif- 
ferential equations which must be solved numerically. 
The first observation is that Eqs.(^ can be viewed as a 
set of Hamilton equations for the complex degrees of free- 
dom. Our approximate dynamics preserves the number 
of particles: N = a^ak, as well as the total energy of 
the system. It therefore corresponds to a genuine micro- 
canonical description. Moreover, the resulting equations 
resemble the famous Fermi-Pasta-Ulam problem of 
a system of harmonic oscillators coupled by a nonlinear 
interaction. A one-dimensional version of this dynamics 
has been studied recently |^ in the context of the pure 
Bose-Einstein condensate (T = 0) reloaded from a har- 
monic into a rectangular trap. Equations studied here, 
however, in spite of a formal analogy, describe quite a 
different physical situation. Our complex amplitudes are 
not expansion coefficients of the condensate wave func- 
tion in some convenient basis. They represent a number 
of coupled "mean fields" - a natural extension of the 
condensate mean field of the Bogolubov approach. Let 
us notice that if we start with 100% occupation of the 
k = mode and keep only this mode in the model we 
simply recover the standard Gross-Pitaevskii equation for 
the interacting condensate. 

For the initial conditions we assume "Bose-Einstein- 
like" occupation of different trap states. While calculat- 
ing the energy in order to determine this distribution, 
however, we neglect the energy of interparticlc interac- 
tions. This is evidently not the equilibrium distribu- 
tion for the interacting system. Moreover, a population 
of individual states does not specify initial conditions 
uniquely. It defines the modulus of the corresponding 
amplitude but says nothing about its phase. In our ap- 
proach each mode is assigned an initial, randomly chosen, 



phase. Any subsequent dynamics depends on the initial 
phases and, in a sense, a single simulation describes a sin- 
gle experimental realization. Microcanonical expectation 
values would require, therefore, an average over these ini- 
tial phases. As we have checked, instead of doing this, 
it is enough to start with some random phases and trace 
the system evolution for a sufficiently long time. The ob- 
served self-averaging can be attributed to the ergodicity 
of the studied dynamics. Contrary to the 1-D analogue 
which is completely integrable , the 3-D version of the 
dynamics may be chaotic. We make use of this fact in 
our simulations and avoid averaging over phases of the 
complex amplitudes. 

Values of the parameters in the model are ^ =71.373 
Hz, N = 10^ and g =0.018 Hz (the atomic mass and the 
scattering length are those of ^^Rb and the size of the 
box is equal to the Thomas- Fermi radius of a condensate 
of N atoms in a trap with frequency of loq = 27r 80 Hz). 
We performed our calculations for the model with 729 
modes (n^ = —4,..., 4, i = x,y,z). Further increasing 
of the number of modes does not lead to a substantial 
change in the results for the case studied in this paper. 
Our calculations show that after a time of the order of 
a few miliseconds the system reaches a dynamical equi- 
librium. The mean occupation of the condensate (k=0 
mode) stabilizes at some value and on larger time scales 
(of the order of a second) it only fluctuates around this 
mean value - see Figure |l|. 
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FIG. 1. Condensate occupation as a function of time for 
total energy per particle E/h = 539Hz. 

The similarity to the Fermi-Pasta-Ulam problem may 
cast some doubts on the genuine ergodicity of the dy- 
namics investigated in this paper . Originally, Fermi, 
Pasta, and Ulam intended to test the ergodic hypothesis 
in the chain of coupled harmonic oscillators. However, 
they observed a quasiperiodic behavior identified by the 
returns of energy to the initial (lowest) mode. This kind 
of behavior has been discovered also by J.H. Eberly and 
co-workers [2^ ] in a resonantly coupled system composed 
of a two-level atom and a single mode of a monochro- 
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matic electromagnetic field being initially in a coherent 
state. Occurrence of the so called revivals in the system 
proves the quantum nature of the electromagnetic field. 
In our calculations the largest time scale for which we 
have studied the dynamics was of the order of one sec- 
ond. On this time scale we did not observe any revivals 
but this, obviously, does not exclude the possibility of re- 
vivals on much larger time scales. In fact our numerical 
simulations involve a finite number of modes and so the 
numerical implementation inevitably leads to a quasiperi- 
odic evolution. 

Since the evolution (see Eqs.(|^)) is Haniiltonian, the 
most natural choice of independent variables are the en- 
ergy and the particle number. The number of particles 
was fixed through all our calculations {N = 10^) and the 
energy of the system was the control parameter. Tra- 
ditionally, however, the temperature, not the energy, is 
used as an independent thermodynamic variable. Calcu- 
lation of the microcanonical temperature requires moni- 
toring of the entropy of the system for different energies. 
Although, in principle, this can be done [p9| , we think 
that the energy is a much better characteristic of the 
system. The reason is twofold. First of all, it is not obvi- 
ous whether the ergodic hypothesis can be applied in the 
studied case and therefore whether the notion of the tem- 
perature can be unambiguously introduced. Secondly, in 
realistic experiments it is rather the final energy of ther- 
mal atoms which is detected in destructive time-of-flight 
measurements. The temperature is a parameter fitted to 
the observed velocity distribution. 
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FIG. 2. Probability distribution of the condensate (k = 
mode) population. Different colors signify different total en- 
ergies per particle: blue E/h = 2036 Hz, green E/h = 1680 
Hz, orange E/h = 1424 Hz, yellow E/h = 1164 Hz, red 
E/h = 539 Hz. 

In Figure || we present probability distributions of the 
condensate population. The different colors signify dif- 
ferent total energies. Let us stress two main features 
which can be seen in Figure ^j: (i) positions of the max- 
ima of these distributions correspond to a mean (time- 
averaged) occupation of the Bose-Einstein condensate. 



(ii) the spread of these distributions around the most 
probable value is a measure of condensate fluctuations. 
The figure clearly indicates that the condensate popula- 
tion decreases monotonically with energy while fluctua- 
tions become larger. 

The occupation and the fluctuations of the interacting 
condensate are shown in detail in Figure ^, where the 
mean occupation of the condensate is depicted (in blue). 
We see that the condensate disappears at the energy per 
particle close to E/h = 2000 Hz. 
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FIG. 3. Condensate occupation (in blue) and fluctuations 
(in red) versus total energy per particle. 

Figure || also presents fluctuations of the condensate 
(in red). Both the mean occupation and the fluctuations 
of the condensate are smooth functions of the energy. 
They do not show any discontinuity signifying a phase 
transition because they correspond to a finite system of 
N = 10^ atoms. Let us notice, however, that fluctua- 
tions reach the maximum value at the energy per par- 
ticle close to E/h — 1450iJz. Moreover, this value of 
energy corresponds to the inflection point of the mean 
occupation of the condensate. Both curves distinguish 
the same characteristic value of the energy. Close to this 
energy the system undergoes rapid changes. This charac- 
teristic energy corresponds to the critical energy for the 
Bose-Einstein condensation. 

In conclusions: flrst of all, we have proposed a very 
efficient approach to the time evolution of an interacting 
Bose-Einstein condensate. This approach is valid for a 
wide range of temperatures except for very small ones. 
We have shown that the method works very well in a 
realistic and relevant range of parameters: energy, num- 
ber of particles, and interaction strength. We calculated 
the mean occupation as well as the fiuctuations of the 
interacting condensate by averaging the corresponding 
time-dependent quantities. We did not explore all pos- 
sible applications of our method. Our aim was rather 
to show its potential. We believe that the semi-classical 
approach presented here is perfectly suited for studying 
many properties of an interacting condensate which so 
far, due to their complexity, were beyond direct theoret- 
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ical investigations. In particular, it may shed new light 
at the unresolved problem of nucleation of a condensate 
and other dynamical phenomena. The method is easily 
extendible to more than one field which allows to study 
coupled (via photoassociation or Feshbach resonances) 
atomic and molecular systems (see |30 ) or other mix- 



tures of Bose gases at finite temperatures. 
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